﻿#include "MeshModifier.h"
#include <vtkm/Bounds.h>
#include <vtkm/Geometry.h>
#include <vtkm/TypeTraits.h>
#include <vtkm/worklet/WorkletCellNeighborhood.h>
#include <vtkm/worklet/WorkletMapTopology.h>

#include <vtkm/Matrix.h>
#include <vtkm/exec/ParametricCoordinates.h>
#include <vtkm/worklet/WorkletMapTopology.h>
#include "Types.h"
#include "Mesh.h"

namespace Loops
{
struct MeshModifier_ForwardStep : vtkm::worklet::WorkletCellNeighborhood
{
  using ControlSignature = void(CellSetIn,
                                FieldInNeighborhood centroid,
                                FieldOut blk_id,
                                FieldOut bnd_id);
  using ExecutionSignature = void(Boundary, _2, _3, _4);

  template<typename NeighIn, typename BndIdType>
  VTKM_EXEC void operator()(const vtkm::exec::BoundaryState& boundary,
                            const NeighIn& centroid,
                            vtkm::Int16& blk_id,
                            BndIdType& bnd_id) const
  {
    blk_id = 0;
    if (!IsCellValid(centroid.Get(0, 0, 0)))
    {
      blk_id = -1; // 不可用单元
    }

    bnd_id = { -1, -1, -1, -1 };
    if (boundary.IJK[0] == 0) // left, inlet
    {
      bnd_id[3] = 3;
    }
    if (boundary.IJK[0] == boundary.PointDimensions[0] - 1) // right, outlet
    {
      bnd_id[1] = 1;
    }
    if (boundary.IJK[1] == boundary.PointDimensions[1] - 1) // top, wall
    {
      bnd_id[2] = 2;
    }
    if (boundary.IJK[1] == 0) // bottom
    {
      bnd_id[0] = 0;
    }

    if (!IsCellValid(centroid.Get(-1, 0, 0)))
    {
      bnd_id[3] = 4;
    }
    if (!IsCellValid(centroid.Get(1, 0, 0)))
    {
      bnd_id[1] = 4;
    }
    if (!IsCellValid(centroid.Get(0, -1, 0)))
    {
      bnd_id[0] = 4;
    }
    if (!IsCellValid(centroid.Get(0, 1, 0)))
    {
      bnd_id[2] = 4;
    }
  }

  VTKM_EXEC bool IsCellValid(Vec3f p) const
  {
    //return Square(p);
    return Sphere(p);
    //return Naca0012(p);
    //return ForwardStep(p);
  }

  VTKM_EXEC bool Square(const Vec3f& p) const
  {
    Real x = p[0];
    Real y = p[1];
    Real a = 0.1;
    if (x < a && x > -a && y < a && y > -a)
    {
      return false;
    }
    else if (x < 0.4 && x > 0.3 && y < 0.7 && y > 0.6)
    {
      return false;
    }
    else if (x < 0.8 && x > 0.6 && y < 0.4 && y > 0.2)
    {
      return false;
    }
    else if (x < 1.5 && x > 1.3 && y < 0.2 && y > -0.1)
    {
      return false;
    }
    else
    {
      return true;
    }
  }

  VTKM_EXEC bool Sphere(const Vec3f& p) const
  {
    if (vtkm::Magnitude(p) < 0.1)
    {
      return false;
    }
    /*if (vtkm::Magnitude(p - Vec3f{ 0.5, 0.2, 0 }) < 0.2)
    {
      return false;
    }
    if (vtkm::Magnitude(p - Vec3f{ 0.5, -0.2, 0 }) < 0.2)
    {
      return false;
    }
    if (vtkm::Magnitude(p - Vec3f{ 0.7, 0.2, 0 }) < 0.15)
    {
      return false;
    }
    if (vtkm::Magnitude(p - Vec3f{ 0.7, -0.2, 0 }) < 0.15)
    {
      return false;
    }
    if (vtkm::Magnitude(p - Vec3f{ 0.9, 0.3, 0 }) < 0.15)
    {
      return false;
    }
    if (vtkm::Magnitude(p - Vec3f{ 0.9, -0.3, 0 }) < 0.15)
    {
      return false;
    }
    if (vtkm::Magnitude(p - Vec3f{ 1.7, -0.3, 0 }) < 0.15)
    {
      return false;
    }
    if (vtkm::Magnitude(p - Vec3f{ 1.5, 0.5, 0 }) < 0.15)
    {
      return false;
    }
    if (vtkm::Magnitude(p - Vec3f{ 2.5, 0.3, 0 }) < 0.15)
    {
      return false;
    }
    if (vtkm::Magnitude(p - Vec3f{ 2.5, -0.3, 0 }) < 0.3)
    {
      return false;
    }*/

    return true;
  }

  VTKM_EXEC bool Naca0012(const Vec3f& p) const
  {
    Real x = p[0];
    Real y = p[1];

    Real a0 = +0.2969, a1 = -0.1260, a2 = -0.3516, a3 = +0.2843, a4 = -0.1036;
    Real t = 0.12;

    if (x < 0 || x > 1)
    {
      return true;
    }
    else
    {
      Real yt =
        5 * t * (a0 * vtkm::Sqrt(x) + a1 * x + a2 * x * x + a3 * x * x * x + a4 * x * x * x * x);
      if (y > yt || y < -yt)
      {
        return true;
      }
      else
      {
        return false;
      }
    }
  }

  VTKM_EXEC bool ForwardStep(const Vec3f& p) const
  {
    vtkm::Bounds box(Vec3f{ 0.6, 0, 0 }, Vec3f{ 3, 0.2, 0 });
    if (box.Contains(p))
    {
      return false;
    }
    return true;
  }
};

void MeshModifier(Mesh& mesh)
{
  vtkm::cont::Invoker invoker;
  const auto& centroid = mesh.GetFieldArray<Vec3f>("center");
  ArrayHandle<vtkm::Int16> blk_id;
  ArrayHandle<Vec<vtkm::Int16, 4>> bnd_id;

  vtkm::cont::Invoker{}(MeshModifier_ForwardStep{}, mesh.GetCellSet(), centroid, blk_id, bnd_id);

  mesh.GetDataSet().AddCellField("blk_id", blk_id);
  mesh.GetDataSet().AddCellField("bnd_id", bnd_id);
}
}